#ifndef __SPLINE_TEST_HPP__
#define __SPLINE_TEST_HPP__

#include "Spline.hpp"
#include <fstream>
using namespace std;

void plot_Xbasis(vector<double> _X, string filename){
    ofstream outfile;
    outfile.open(filename);
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << _X[i] << ",";
    outfile << "]\n";
    outfile.close();
}

void plot_dots(vector<double> _X, string filename){
    ofstream outfile;
    outfile.open(filename, ios::app);
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << _X[i] << ",";
    outfile << "]\n";
    outfile.close();
}

void plot_Func(func& F, vector<double> _X, string filename){
    ofstream outfile;
    outfile.open(filename, ios::app);
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << F(_X[i]) << ",";
    outfile << "]\n";
    outfile.close();
}

void plot_spline(SplineBase& spline, vector<double> _X, string filename){
    ofstream outfile;
    outfile.open(filename, ios::app);
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << spline(_X[i]) << ",";
    outfile << "]\n";
    outfile.close();
}

void plot_curve(CubicSpline_2dcurve& spline, vector<double> _X, string filename){
    ofstream outfile;
    outfile.open(filename);
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << spline(_X[i]).x() << ",";
    outfile << "]\n";
    outfile << "[";
    for(int i=0;i<_X.size();i++) outfile << spline(_X[i]).y() << ",";
    outfile << "]\n";
    outfile.close();
}

class f_in_A:public func{
public:
    double operator () (double x){
        return 1/(1+25*x*x);
    }
    double diff(double x){
        return -50*x/(1+25*x*x)/(1+25*x*x);
    }
};

class f_in_B:public func{
private:
    double shift;
public:
    f_in_B(double shift):shift(shift){}
    ~f_in_B(){}
    double operator()(double y){
        double x = y+shift;
        return 1/(1+x*x);
    }
    double diff(double y){
        double x = y+shift;
        return -2*x/(1+x*x)/(1+x*x);
    }
};

class _pow:public func{
private:
    int _pow_=1;
public:
    _pow(int pow):_pow_(pow){}
    ~_pow(){}
    double operator() (double x){
        return pow(x,_pow_)/pow(6,_pow_)*0.128924;
    }
};

void TaskA(string filename){
    int N[] = {6, 11, 21, 41, 81};
    vector<double> _N;
    for(int i=0;i<5;i++) _N.push_back(N[i]);
    vector<double> X;
    vector<double> _maxnorm;
    int _Div = 1000;
    double step = 2.0/_Div; double x = -1;
    for(int i=0;i<=_Div;i++){
        X.push_back(x);
        x += step;
    }
    //plot_Xbasis(X, filename+".dat");
    f_in_A F;
    for (int i=0;i<5;i++){
        plot_Xbasis(X, filename+"_"+to_string(N[i])+".dat");
        step = 2.0/N[i]; x = -1;
        vector<double> knots;
        for(int j=0;j<=N[i];j++){
            knots.push_back(x);
            x += step;
        }
        //vector_print(knots);
        CubicSpline_func Spline(knots, F, 2);
        double M = 0;
        for(int i=0;i<knots.size()-1;i++){
            double _mid = (knots[i]+knots[i+1])/2;
            double temp = abs(Spline(_mid)-F(_mid));
            if(M < temp) M = temp; 
        }
        cout << "The maxnorm is " << M << " when N equals " << N[i] << "." << endl;
        _maxnorm.push_back(M);
        plot_Func(F, X, filename+"_"+to_string(N[i])+".dat");
        plot_spline(Spline, X, filename+"_"+to_string(N[i])+".dat");
    }
    vector<double> _logx;
    for(int i=5;i<=100;i++) _logx.push_back(i);
    plot_Xbasis(_logx,filename+"_Log.dat");
    for(int i=2;i<=6;i++){
        _pow _Po(-i);
        plot_Func(_Po, _logx,filename+"_Log.dat");
    }
    plot_dots(_N,filename+"_Log.dat");
    plot_dots(_maxnorm,filename+"_Log.dat");
}

void TaskC(string filename){
    vector<double> t1;
    vector<double> t2;
    int N1 = 11; int N2 = 10;
    for(int i=1;i<=N1;i++) t1.push_back(i-6);
    for(int i=1;i<=N2;i++) t2.push_back(i-5.5);

    vector<double> X1; vector<double> X2;
    vector<double> _X1; vector<double> _X2;
    int _Div = 1000; 
    double step1 = (t1[t1.size()-1]-t1[0])/_Div;
    double step2 = (t2[t2.size()-1]-t2[0])/_Div;
    double x1 = t1[0]; double x2 = t2[0];
    for(int i=0;i<=_Div;i++){
        X1.push_back(x1); _X1.push_back(x1+6); x1 += step1;
        X2.push_back(x2); _X2.push_back(x2+5.5); x2 += step2;
    }

    f_in_B F(-6);
    BSpline_cubic_func BS1(N1, F);
    BSpline_quad_func BS2(N2, F);

    plot_Xbasis(X1,filename+"_cubic.dat");
    plot_Func(F, _X1, filename+"_cubic.dat");
    plot_spline(BS1, _X1, filename+"_cubic.dat");

    plot_Xbasis(X2,filename+"_quad.dat");
    plot_Func(F, _X2, filename+"_quad.dat");
    plot_spline(BS2, _X2, filename+"_quad.dat");
}

void TaskD(){
    vector<double> t1;
    vector<double> t2;
    int N1 = 11; int N2 = 10;
    for(int i=1;i<=N1;i++) t1.push_back(i-6);
    for(int i=1;i<=N2;i++) t2.push_back(i-5.5);

    vector<double> X1; vector<double> X2;
    vector<double> _X1; vector<double> _X2;
    int _Div = 1000; 
    double step1 = (t1[t1.size()-1]-t1[0])/_Div;
    double step2 = (t2[t2.size()-1]-t2[0])/_Div;
    double x1 = t1[0]; double x2 = t2[0];
    for(int i=0;i<=_Div;i++){
        X1.push_back(x1); _X1.push_back(x1+6); x1 += step1;
        X2.push_back(x2); _X2.push_back(x2+5.5); x2 += step2;
    }

    f_in_B F(-6);
    BSpline_cubic_func BS1(N1, F);
    BSpline_quad_func BS2(N2, F);

    double _cal[] = {-3.5,-3,-0.5,0,0.5,3,3.5};
    double E1 = 0; double E2 = 0;
    cout << "---- x ----|----- quad_err -----|----- cubic_err -----|" << endl;
    for(int i=0;i<7;i++){
        double x = _cal[i];
        double err1 = abs(BS2(x+6)-F(x+6));
        double err2 = abs(BS1(x+6)-F(x+6)); 
        printf("%+11f|%+20f|%+21f|\n",x,err1,err2);
        E1 += err1; E2 += err2;
    }
    cout << "The sums of error are respectively " << E1 << " and " << E2 << "." << endl;
}

void TaskE(string filename){
    int _N[]={10,40,160};
    double epsilon = 0.000001;
    double a = -sqrt(3)+epsilon; double b = sqrt(3)-epsilon;
    int _Div = 100;
    vector<double> Xbasis;
    double x = 0;
    for(int i=0;i<=_Div;i++){
        Xbasis.push_back(x);
        x += 1.0/_Div; 
    }
    for(int i=0;i<3;i++){
        vector<double> X;
        vector<double> Y;
        int N = _N[i]/2;
        double step = (b-a)/N;
        x = a;
        for(int j=0;j<=N;j++){
            X.push_back(x);
            Y.push_back((-sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x += step;
        }
        x = b;
        for(int j=0;j<=N;j++){
            X.push_back(x);
            Y.push_back((sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x -= step;
        }
        CubicSpline_2dcurve Curve(X,Y);
        //vector_print(Xbasis);
        //vector_print(X);
        //vector_print(Y);
        plot_curve(Curve, Xbasis, filename+"_"+to_string(_N[i])+".dat");
    }
}

void TaskEOPT(string filename){
    int _N[]={10,40,160};
    double epsilon = 0.000001;
    double a = -sqrt(3)+epsilon; double b = sqrt(3)-epsilon;
    int _Div = 100;
    vector<double> Xbasis;
    double x = 0;
    for(int i=0;i<=_Div;i++){
        Xbasis.push_back(x);
        x += 1.0/_Div; 
    }
    for(int i=0;i<3;i++){
        vector<double> X;
        vector<double> Y;
        int N = _N[i]/2;
        double step = (b-0)/(N/2);
        x = 0;
        for(int j=0;j<N/2;j++){
            X.push_back(x);
            Y.push_back((-sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x += step;
        }
        step = (b-0)/(N-(N/2));
        for(int j=0;j<N-N/2;j++){
            X.push_back(x);
            Y.push_back((sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x -= step;
        }
        step = (0-a)/(N-(N/2));
        for(int j=0;j<N-N/2;j++){
            X.push_back(x);
            Y.push_back((sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x -= step;
        }
        step = (0-a)/(N/2);
        for(int j=0;j<N/2;j++){
            X.push_back(x);
            Y.push_back((-sqrt(3-x*x)+sqrt(abs(x)))*2/3);
            x += step;
        }
        X.push_back(0);
        Y.push_back(-1.1547);
        CubicSpline_2dcurve Curve(X,Y);
        //vector_print(Xbasis);
        //vector_print(X);
        //vector_print(Y);
        plot_curve(Curve, Xbasis, filename+"_"+to_string(_N[i])+".dat");
    }
}

void TestLinear(string filename){
    vector<double> t1;
    vector<double> t2;
    int N1 = 11; int N2 = 10;
    for(int i=1;i<=N1;i++) t1.push_back(i-6);
    for(int i=1;i<=N2;i++) t2.push_back(i-5.5);

    vector<double> X1; vector<double> X2;
    int _Div = 1000; 
    double step1 = (t1[t1.size()-1]-t1[0])/_Div;
    double step2 = (t2[t2.size()-1]-t2[0])/_Div;
    double x1 = t1[0]; double x2 = t2[0];
    for(int i=0;i<=_Div;i++){
        X1.push_back(x1); x1 += step1;
        X2.push_back(x2); x2 += step2;
    }
    f_in_B F(0);
    vector<double> f1;
    vector<double> f2;
    for(int i=0;i<t1.size();i++) f1.push_back(F(t1[i]));
    for(int i=0;i<t2.size();i++) f2.push_back(F(t2[i]));
    Linear_Spline L1(t1,f1);
    Linear_Spline L2(t2,f2);

    plot_Xbasis(X1,filename+"_1.dat");
    plot_Func(F, X1, filename+"_1.dat");
    plot_spline(L1, X1, filename+"_1.dat");

    plot_Xbasis(X2,filename+"_2.dat");
    plot_Func(F, X2, filename+"_2.dat");
    plot_spline(L2, X2, filename+"_2.dat");

}

#endif